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Abstract 

This tutorial presents an introduction into continuum descriptions of cytoskeletal dynamics. In contrast to discrete 
models in which each molecule keeps its identity, such descriptions are given in terms of averaged quantities per 
unit volume like the number density of a certain molecule. Starting with a discrete description for the assembly 
dynamics of cytoskeletal filaments, we derive the continuity equation, which serves as the basis of many 
continuum theories. We illustrate the use of this approach with an investigation of spontaneous cytoskeletal 
polymerization waves. Such waves have by now been observed in various cell types and might help to orchestrate 
cytoskeletal dynamics during cell spreading and locomotion. Our analysis shows how processes at the scale of 
single molecules, namely, the nucleation of new filaments and filament treadmilling, can lead to the spontaneous 
appearance of coherent traveling waves on scales spanning many filament lengths. For readers less familiar with 
calculus, we include an informal introduction to the Taylor expansion. 



Introduction 

With the advancements in microscopy techniques it has 
become increasingly clear that a true understanding of 
many cellular phenomena requires to take spatial aspects 
into account. This is clearly the case for the cytoskeleton 
as illustrated by the changes in the microtubule network 
during cell division or the reorganization of the actin 
meshwork during cell locomotion. To reach a quantita- 
tive understanding of the mechanisms underlying these 
processes, concepts and methods from physics can be 
extremely valuable. These methods include notably the 
theoretical analysis of cellular systems. Most biological 
and medical curricula today lack, unfortunately, a thor- 
ough introduction into mathematical and physical tools, 
which often leads to discomfort on the part of biologists 
and physicians, when confronted with the results of a 
theoretical study. This holds in particular for continuum 
descriptions. This tutorial is meant to familiarize life 
scientists with the basic ideas underlying this approach. 

Continuum theories have had a great success in 
describing static and dynamic phenomena for large 
classes of matter. Well-known examples range from 
simple fluids [1] and elastic materials [2] to fluid 
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membranes and vesicles [3]. Less conventional examples 
include the cytoskeleton [4] and even flocks of birds [5] . 
Especially for flocks of birds, one might at first sight be 
inclined to use a discrete approach instead, in which the 
position and behavior of each individual bird is consid- 
ered. Such models indeed exist and have yielded valu- 
able insights, see for example [6]. Similarly, simulation 
tools like cytosim allow the user to study cytoskeletal 
dynamics, while keeping track of each individual cytos- 
keletal filament, of each motor molecule, and of any 
other cytoskeletal protein possibly present, for example, 
the passive cross-linker a-actinin [7]. There are, how- 
ever, two major problems associated with this approach. 
One is most easily revealed by considering water. In a 
discrete description, we would need to track the position 
of each individual water molecule, its momentum, and 
its interaction with all surrounding water molecules. 
Given the astronomic number of water molecules even 
in small volumes, this is a daunting task. Furthermore, 
one needs to provide details of the interactions between 
the individual entities and it is a priori not clear, which 
of these details will be important for the large scale 
dynamics of a system. 

In contrast, in a continuum theory, the ultimately par- 
ticulate nature of a material or system is neglected. 
Intuitively, this approximation is valid for structures on 
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length scales that are very much larger than the discrete 
constituents of the material under investigation. Maybe 
somewhat surprisingly, though, they often describe a 
material's behavior correctly even down to length scales 
that approach the constituents' sizes. The Navier- Stokes 
equation, for example, describes the flow of water in 
micro- and to some extent even in nanofluidic devices 
[8]. From a more technical point of view, continuum 
theories often have the advantage over more micro- 
scopic descriptions to provide a comprehensive account 
of possible material behavior. Furthermore, many of the 
details of the interactions between the constituents are 
lumped together into a comparatively small set of para- 
meters. In this tutorial, we show how continuum the- 
ories can be developed and used in a cellular context to 
describe cytoskeletal dynamics. 

Our aim is not to enable you to formulate and solve 
your own continuum theories. Still, you should get a 
feeling for the method and for what you can and cannot 
expect from this approach. This tutorial necessarily con- 
tains a certain amount of mathematics. To make the 
tutorial accessible to an audience not using calculus on 
an everyday basis, we have tried to make all steps in the 
calculations explicit. You should have a pencil and a 
sheet of paper ready to follow them in detail. In addi- 
tion, some calculations are left for the reader. To per- 
form them you do not need much beyond a sincere 
interest in these techniques and some basic knowledge 
of calculus, notably to take derivatives. To keep the 
mathematics tractable, we have also simplified certain 
biological processes. The emphasis of this tutorial is 
really on mathematical and physical concepts rather 
than on a detailed analysis of cellular processes. 

The tutorial is organized as follows. We start with a 
discrete description of the length dynamics of cytoskele- 
tal filaments from which we derive a continuum theory. 
This approach will lead us to the continuity equation. 
We will use this equation to develop a continuum 
description of treadmilling filaments interacting with 
nucleation promoting factors. We will use this formal- 
ism to study a possible mechanism underlying sponta- 
neous polymerization waves. Two excursions aim at 
recalling the Taylor expansion and give a general 
account of the continuity equation. 

Mathematical prelude: The Taylor expansion 

A mathematical tool that we will make frequent use of 
when going to the continuum limit is the so-called Tay- 
lor expansion of a (real-valued) function. The aim here 
is not to give a mathematically rigorous account of the 
Taylor expansion, but rather to introduce the basic idea. 
Consequently, we assume that the functions we will deal 
with are all sufficiently well-behaved such that all opera- 
tions we perform on them are allowed. 



The aim of the Taylor expansion is to approximate a 
function in the vicinity of a chosen point by a polyno- 
mial, usually of low order. Consider the function / dis- 
played on Figure la and the point Xq. The value of the 
function at Xq is denoted by f{xo). A very crude approxi- 
mation of /in the vicinity of IS J{xo). An obviously bet- 
ter approximation of /is given by the tangent to /in Xq. 
If we denote the slope of /in xq by f{xo), then the tan- 
gent is given by the expression /^Vg) + /(^o)(* - ^o)> where 
the slope of /in Xq is just the derivative of / in Xq, that is. 



Here, we have suppressed the dependence on Xq for 
the ease of notation. 

Now, the derivative of a function is itself a function. If 
we use the tangent to approximate the function / then 
this amounts to the same as to approximating the deri- 
vative by a constant, namely /(xq), see Figure lb. As 
above, we can improve the approximation by using the 
tangent to / at Xq- Denoting the slope of the derivative 
in Xo by /'(xq), we can write f{x) ~f(xo) + f'{xo)(x - Xq)- 
To be consistent, we should also modify our approxima- 
tion for / Putting 

/(x) ^ f[Xo) + /'(xo)(x - Xo) + ^/"(xo)(x - Xo)^, (2) 

we see that the derivative of the approximation for/ 
equals the approximation for /. The degree of the poly- 
nomial used to approximate a function /is called the 
order of the approximation. 

Your turn: Calculate the Taylor expansion of the 
exponential function expix) with respect to Xq = 0 up to 
second order. 

Obviously, this scheme can be iterated further to get 
better and better approximations of /in the vicinity of 

Xq. However, we will content ourselves in the following 
with approximations up to second order at most. 

Length-dynamics of active polar filaments - part I 

As a first example for a continuum description of cytos- 
keletal dynamics, we will consider the assembly and dis- 
assembly of cytoskeletal filaments. Actin filaments and 
microtubules are linear aggregates of non-covalently 
bound subunits, namely G-actin and tubulin dimers, 
respectively. G-actin is a molecule of about 42 kDa with 
linear extensions of about 6.7 nm and 4 nm, whereas a 
tubulin dimer has a molecular weight of about 110 kDa 
and is about 8 nm long and 4 nm in diameter. The 
length distributions of F-actin and of microtubules in 
cells is currently largely unknown, but typical values 
have been estimated to be around 500 nm up to a few 
microns for F-actin and up to several tens of microns 
for microtubules. These subunits are structurally polar 
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Figure 1 The Taylor expansion, a) Various polynomial 
approximations to a function fin the vicinity of x = Xq (blue line). The 
simplest approximation is to use the value f(xo) (dotted line). The 
"first-order approximation" is given by the tangent to f in x = Xq (green 
ine). The "second-order approximation" is given by the quadratic 

polynomial /(Xq) +f'{Xo){x — Xq) + jf"{XQ){x — Xq)^ 

(red line), b) Polynomial approximations to the derivative f of fin 
the vicinity of x = Xq (green line). The tangent to f in x = Xq is 
given by the derivative of the second-order approximation to 
f (red line). 



resulting in polar filaments with different subunit 
attachment and removal rates at the two ends. For 
further reference, we will call the fast growing end the 
plus-end and the other end the minus-end. In addition, 
G-actin and tubulin can bind nucleotides and their rates 
of attachment to and removal from a filament depend 
on the state of the nucleotide bound, that is, whether 
it is a nucleoside triphosphate (NTP) or a nucleoside 
diphosphate (NDP). Through hydrolysis of an NTP 
energy can be introduced into the system on a mole- 
cular level. This is the defining property of an active 
material. 

The structural polarity together with the activity 
makes the assembly and disassembly dynamics of 



cytoskeletal filaments more interesting than that of most 
polymers that are usually studied. In particular, it can 
lead to a phenomenon called treadmilling, where the 
plus-end grows on average, while the minus-end shrinks 
on average [9,10]. This is the situation we will consider 
in the following without further questioning the origin 
of this phenomenon. Also, we will consider generic 
aspects of treadmilling filaments rather than specifically 
actin filaments or microtubules, because here we are 
mostly interested in introducing the concept of conti- 
nuum theories. 

We consider the situation of active polar filaments 
that assemble in a solution containing subunits. 
Let Ci denote the concentration of filaments consisting 
of i = 1, 2, ... subunits. For now, we assume that the 
filaments are homogeneously distributed and will con- 
sider a possible space-dependence of the distribution 
later. The time evolution of Cl is given by 

-^Cl = feflCi-i - feaCi - IZdCi + kdCL+l, (3) 

at 

where ka and are the rates of subunit attachment at 
the plus- and of subunit removal at the minus-end, 
respectively. These rates depend strongly on the envir- 
onmental conditions and on the state of the nucleoside 
bound, but typical values are about lOs'^^fiM'^ for the 
attachment rate of ATP-G-actin at the plus-end of an 
actin filament and 0.25s''^ for the removal rate of ADP- 
G-actin at the minus-end [11]. Note, however, that the 
effective rates of subunit addition at the plus- and subu- 
nit removal at the minus-end result from complex 
dynamic processes inside the filament, such that they 
can strongly deviate from these values. For the time 
being, we will consider them to be constant, which 
amounts to assuming a constant reservoir of filament 
subunits. We will discuss a possible experimental reali- 
zation below. 

While there is in principle nothing wrong with this 
discrete description, it turns out that a continuum 
description is easier to analyze. The idea is to neglect 
the discrete nature of the subunits and consider a con- 
tinuous length distribution. Indeed, if you are interested 
in the behavior of the length distribution on scales that 
are larger than many subunit sizes, many details are 
irrelevant. 

Technically, we start by replacing the concentrations 
Cl by a density c{t), such that = c{LS)S, where S is the 
change in filament length upon addition or removal of a 
subunit. We thus have 

^c(L«) = fe„c((L - 1)5) - fe„c(LS) - }idc(a) + kjc[{L + 1)S). (4) 

If we are interested in structures that are significantly 
larger than d, we can treat (5 as a small quantity. 
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Consequently, we can use the Taylor expansion to 
approximate c({L ± 1)S) by a linear function: 

c{[L±l)S ^ c{LS)±c'{LS)S. (5) 

Inserting this expression in Eq. (4) we obtain 

Yt:{LS) = h\c{LS] - c'{LS]S] - k„c{LS) - hc{LS) + ki\c{LS) + cf{LS]S] (6) 



-kad{LS)S + kdd{LS)S 



(7) 



The important step is now to replace the discrete 
expression Ld by the continuous variable i. Intuitively, 
this makes sense as long as you do not consider features 
of the distribution c(€) on length scales that are of the 
order of the subunit size or smaller. It also implies that 
molecular details of the assembly processes can be irre- 
levant for the large scale structures. We thus arrive at 



3 9 

-C(£,t) = -{Va-Vj)—C{l,t). 



(8) 



Here, = kaS and = k^S are the assembly and dis- 
assembly velocities. Furthermore, in this expression the 
derivative with respect to time, d/dt has been replaced 
by a partial derivative d/dt. The reason is that the den- 
sity c depends on two continuous variables, £ and t, and 
taking the partial derivative indicates that we want to 
calculate the derivative with respect to t for fixed £. 
Analogously, d/di indicates that we take the derivative 
with respect to € for fixed t. 

Interlude: The continuity equation 

The partial differential equation that we derived for the 
length dynamics of treadmilling filaments can be cast 
into the following form 



9 . . 9 



(9) 



where j{£, t) = [v^ - v^)c{£, t). Such an equation is 
known as a continuity equation. It generally describes 
the evolution of a conserved quantity with / being the 
associated current. We will now see, how such an equa- 
tion can be obtained on general grounds. 

For simplicity, we start with a one-dimensional system 
and will then consider the general case. Consider a 
quantity Q that is distributed along an axis. This could, 
for example, be the electric charge along a DNA 
molecule. Let's denote by Q{x) the amount of Q in the 
intervall {x, x + Ax] and assume that Q can neither be 
generated nor destroyed. In that case Q(x) can only 
change through a flux of Q across the boundaries of the 
interval. Let's denote this flux by which by convention 
is positive if directed towards larger x, see Figure 2a. 
Formally, we then have the flux balance equation 



a) 



j(x+dx) 



X 



x+dx 



j(x,y+dy) 



j{x+dx,y) 



b) 



y+dy- 



y 



X X I dx 

Figure 2 The continuity equation, a) Flux balance in one 
dimension: the quantity Q can enter or leave the interval [x, x + /\x] 
only via flux j through the boundaries, b) Flux balance in two 
dimensions: the quantity Q can enter or leave the rectangular 
region of size Ax x Ay only via flux j through the boundaries. 




— Q(x) =i{x) — i(x+ Ax), 
at 

Using the same expansion as above, we find 



d 

-— j(x)Ax. 
ax 



3 



(10) 



(11) 



(12) 



If we divide by Ax and consider the limit Ax — > 0, 
then we get 



3t 



3 

dx 



(13) 



where q(x) = limA^c^o Q(x)/Ax is the (line) density 
associated with Q. This is the one-dimensional continu- 
ity equation and has the same form as Eq. (9). 

If we go to higher dimensions, then the flux becomes 
a vector j with components jx> jy> and so on denoting 
the components of the current in x- and j-direction, 
respectively. To calculate the change of a quantity in a 
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small box, one then has to consider the contributions of 
the current in all directions, see Figure 2b. In two spatial 
dimensions, this leads to 



9 _ 9 9 



(14) 



Note, that here q is now a surface density. 

Your turn: Derive Equation (14). How does it look in 
three spatial dimensions? 

The right hand side of Eq. (14) is called the divergence 
of j and is often written as V-j. Adding a source and sink 
term Sq to account for possible generation or destruc- 
tion of Q, we can write 



9 



Sq, 



(15) 



which is the most general form of the continuity 
equation. 

Example: A reaction-diffusion system 

As a specific example consider the diffusion of two 
molecular species A and B that can be converted into 
each other. For specificity you might have a protein in 
mind that can change between two different conforma- 
tions. We describe the distribution of the molecules by 
the number densities a and b, respectively. The continu- 
ity equations for the two molecular species then read: 



dt 



(16) 



(17) 



The molecules are assumed to diffuse in an aqueous 
solution, such that their currents are directed along the 
corresponding concentration gradients. Explicitly, the 
components of the currents are (in three dimensions) 



kx = -Da— a 
9x 



^ 9 
-Da — a 

9y 



9 

-Da— t 

9z 



Ja,z ■ 

and 



ih,x = -Db—b 
dx 



(18) 
(19) 
(20) 

(21) 
(22) 



ih,z = -Dh—b. 

az 



(23) 



The diffusion constants and Dj, are positive and 
the minus-signs indicate that the flux is directed from 
high to low concentrations. Since we assume that spe- 
cies A is only generated and destroyed by conversion 
into B and vice versa, the source terms in Eqs. (16) and 
(17) have to have opposite signs. Their functional form 
depends on the molecular details of the conversion pro- 
cesses. If, for example, conversions occur spontaneously 
at a rate from A to B and at rate ojt, from B to A, 
then one would write S = -cOaCi + (Obb- 

Inserting the expressions for the currents and the 
source terms and shifting the terms containing the cur- 
rents to the right hand side, we then get 



— a = DaAa — cOaU + ayi,b 
dt 



dt 



DfcAi? + coaa — (Obb. 



(24) 



(25) 



Here, A is the so-called Laplace operator. In absence 
of the source terms, each of the two equations is the 
common diffusion equation. In presence of the source 
terms, the two equations define a simple reaction-diffu- 
sion system as introduced by A. Turing [12]. 

Your turn: Show that the Laplace operator can be 
written as 



^ _ 92 ^ 92 ^ 92 
9x2 9y2 ■ 



(26) 



Boundary conditions 

To complete the description of the evolution of Q in a 
finite domain one has to specify so-called boundary con- 
ditions. They determine the behavior of the current 
and/or the density at the domain boundaries. For exam- 
ple, if you consider particles in a vessel with impene- 
trable walls, then the current through these walls must 
vanish, jx = 0, where denotes the component of the 
current perpendicular to the wall. 

Length-dynamics of active polar filaments - part II 

Let us return to the length distribution of treadmilling 
filaments and focus on the case that Va > V4 such that 
filaments grow on average. If and are constants, 
then filament growth would be unbounded. In a real 
system, though, the pool of subunits is exhausted at 
some point, which leads to a stop of filament growth. 
We will consider here a formally simpler alternative, 
namely that filaments vanish at a constant rate inde- 
pendent of their length. To our knowledge, there is no 
experimental evidence for such disassembly of neither 
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actin filaments nor microtubules. However, a more rea- 
listic description would at this point only distract from 
our main goal which is to introduce the idea of conti- 
nuum theories. One could, however, realize such a 
dynamics in vitro by removing in regular intervals let's 
say half of the solution containing the filaments and 
subsequently adding new buffer containing filament sub- 
units. We neglect possible breaking of filaments during 
this process as well as filament annealing. 

Since partitioning and distributing filaments results in 
diluting them, we also need to provide a source of new 
filaments. At physiological monomer concentrations, the 
rate at which new actin filaments form spontaneously is 
negligible. Instead cells rely on proteins that assist the 
formation of new filaments, so-called nucleation pro- 
moting factors (NPFs). Common examples of such pro- 
teins are the Arp2/3 complex, which when itself bound 
to an existing filaments generates a new filament that 
has its pointed end connected to the complex, or mem- 
bers of the formin family, which stay attached to the 
barbed end of a newly generated filament and assist 
further elongation. In the following, we will neglect 
these aspects of NPFs and just consider them to gener- 
ate new filaments. 

Let us assume that we maintain a constant density of 
NPFs that generate new filaments at a rate v. We can 
account for this feature by a boundary on the filament 
current at filament length i = 0, that is /[^ = o = The 
length dynamics of polar filaments is then completely 
specified by the following two equations: 



9 9 
9t d£^ " 



{Va - Vd)c{l = 0) = V. 



(27) 



(28) 



From these equations, we can get the distribution of 
filament lengths in steady state, that is, for dc/dt - 0. 
Explicitly, 



Va-Vd 



■ exp{-€/A.} 



(29) 



where X - {va- Vd)lvd is a characteristic length. 

Your turn: Show that the distribution q of Eq. (29) 
indeed solves Eq. (27) with dc/dt = 0 and fulfills the 
boundary condition Eq. (28). 

In the present example, we have derived the conti- 
nuum theory from an underlying discrete microscopic 
picture. In general, different microscopic pictures will 
lead to the same continuum theory. An alternative pos- 
sibility is to start directly with the continuity equation. 
In that case, the definition of the current and the source 
terms reflect essential features of the physics underlying 
the molecular processes of interest. We will give now an 



example of this approach and treat the case of spatially 
heterogeneous distributions of treadmilling filaments 
interacting with NPFs. 

Theoretical description of a network of filaments 

In recent years, actin waves have been observed in var- 
ious organisms [13]. These waves appeared either after 
drug treatment or spontaneously, see Figure 3 for an 
example. It has been suggested that such waves might 
orchestrate the cytoskeletal components for cell migra- 
tion [14,15]. We will now show that the interplay of 
dynamic polar filaments and nucleation promoting fac- 
tors (NPFs) can lead to the spontaneous emergence of 
waves [14,16-18]. The tool we will use in this investiga- 
tion is again a continuum theory. 

To give a full account of the configuration of the actin 
cytoskeleton, we would need to specify the position, 
orientation, length, and configuration of the filaments. 
We will simplify the description by considering only 
length scales larger than the average filament length and 
by treating the filaments as points, while retaining the 
information about their orientation. Furthermore, we 
will consider a situation, where all filaments are aligned 
along a single axis as is the case for actin filaments in 
stress fibers or contractile rings. The processes we 
intend to capture with our description include filament 
treadmilling, filament nucleation, and the dynamics of 
the nucleation promoting factors. They are schematically 
represented in Figure 4. 

Let's turn to the formalism. The state of the filament 
bundle along the ;>c-axis is specified by the densities 
and c of filaments with their plus-end pointing into the 
direction of increasing and decreasing values of x, 
respectively. The time evolution of and c is given by 
continuity equations: 



9 
9£ 



dx'" 



dt 



9 . 

+ /r 



(30) 



(31) 



where and are the respective filament fluxes in 

the ;>(;- direction. The source and sink terms and S' 
account for the generation of new filaments and their 
disappearance. 

We now specify the expressions for the currents and 
the source. As above, we assume treadmilling dynamics 
for the filaments. As a consequence, the filaments will 
move in space into the direction of the filament, see 
Figure 4a. In addition, we will account to some extent 
for fluctuations in the system by introducing an effective 
diffusion current. This current is proportional to gradi- 
ents in the filament density and captures thermal 
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fluctuations and to some extent the stochasticity of 
filament assembly and disassembly. We thus have 

= vc^ - of (32) 
ax 

i~ = —vc~ — D—c~, (33) 

dx 

where v is the treadmilling velocity. As in the previous 
section, we also have to introduce boundary conditions. 
We will use in the following periodic boundary condi- 
tions, that is, for a system of size L, the densities and 
currents are L-periodic functions, for example, j*{x + L) 
= i*{x). We use these boundary conditions mostly for 
computational convenience, but they have also some 
biological relevance, for example, in the case of contrac- 
tile actin rings. Let us recall, however, that in this tutor- 
ial, we do not want to give a detailed account of 
biological processes, but rather focus on formal aspects. 
For a more realistic application of a continuum descrip- 
tion of cytoskeletal waves, see for example Refs. [15-18]. 

The source terms depend on the distribution of nucle- 
ating promoting factors (NPFs). The crucial point is, that 
these have a dynamics on their own, which couples back 
to the distribution of filaments. Motivated by the Arp2/3 
complex, we distinguish two classes of NPFs: bound to 
filaments and cytosolic. Furthermore, we assume that 
bound NPFs generate filaments of the same orientation 
as the filament they are bound to, see Figure 4c. As a 
consequence, we introduce the densities and n of 
NPFs bound to plus- and minus-filaments, respectively, 
and the density ric of cytosolic NPFs. The source terms 
5* and S" then read 

S* = v{rf + enc) - VdC* (34) 
S~ = v[rr + Etic) — VdC~ . (35) 



Here, v is the nucleation rate and e a numerical factor 
that accounts for the weaker nucleation activity of cyto- 
solic NPFs. Finally, similarly to the previous section, 
is the rate of filament disassembly, see Figure 4b. 

Also the dynamics of the NPFs is governed by conti- 
nuity equations. The transport of NPFs is purely diffu- 
sive. This holds also for bound NPFs as they stay 
attached to a specific subunit of a filament and as tread- 
milling does not move filament subunits. Note, however, 
that the system behavior does not change qualitatively if 
NPFs are transported with the filaments, as could be the 
case in presence of molecular motors. For bound NPFs 
we assume a constant unbinding rate m^, see Figure 4e. 
The attachment dynamics of NPFs to filaments is more 
involved. First of all there is a base rate oj^ per filament 
at which cytosolic NPFs bind. In addition, we assume a 
cooperative effect among bound NPFs, namely, that 
bound NPFs favor the binding of additional NPFs, see 
Figure 4d. We refrain from suggesting a possible mole- 
cular picture underlying this effect and simply note that 
some cooperative effect is necessary to spontaneously 
generate waves. Altogether, we arrive at the following 
dynamic equations: 

Also the dynamics of the NPFs is governed by conti- 
nuity equations. The transport of NPFs is purely diffu- 
sive. This holds also for bound NPFs as they stay 
attached to a specific subunit of a filament and as tread- 
milling does not move filament subunits. Note, however, 
that the system behavior does not change qualitatively if 
NPFs are transported with the filaments, as could be the 
case in presence of molecular motors. For bound NPFs 
we assume a constant unbinding rate co^, see Figure 4e. 
The attachment dynamics of NPFs to filaments is more 
involved. First of all there is a base rate cu^ per filament 
at which cytosolic NPFs bind. In addition, we assume a 
cooperative effect among bound NPFs, namely, that 
bound NPFs favor the binding of additional NPFs, see 
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Figure 3 Example of an actin polymerization wave in Dictyostelium discoideum Circular wave expanding the cell border of Dictyostelium 
discoideum. Cells expressing mRFP-LimED to label filamentous actin structures and GFP-Arp3 for incorporation into the Arp2/3 complex were 
recorded by dual-emission TIRF microscopy during recovery from latrunculin A treatment. The protrusion stops when the wave starts to collapse. 
Time is indicated in seconds. Scale ban 10 (jm. From Bretschneider et al. [20] with permission from Biophysical Journal. 
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Figure 4 Schematic representation of filament dynamics in presence of nucleation promoting factors (NPFs) a) Filaments treadmill by 
attachment of monomers at one end and detachment at the other, b) Filaments disassemble spontaneously, c) New filaments are generated by 
NPFs. d) NPFs bind to filaments. Binding is cooperative, favoring attachments in the vicinity of already bound NPFs. e) NPFs spontaneously 
detach from filaments. 



Figure 4d. We refrain from suggesting a possible mole- 
cular picture underlying this effect and simply note that 
some cooperative effect is necessary to spontaneously 
generate waves. Altogether, we arrive at the following 
dynamic equations: 
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at 



1 + coin ) + c (1 + ct>in 



(36) 



(37) 



(38) 



In this expression, is the diffusion constant for 
cytosolic NPFs and coi is a measure of the cooperative 
effect on NPF attachment. We have assumed that NPFs 
are neither generated nor do they degrade. 

Your turn: Show that the total number of NPFs is 
conserved. 

This completes our description of treadmilling 
filaments. There are now several ways to analyze the 
above equations. An analytic solution of such coupled 
non-linear partial differential equations is usually not 
possible and one thus relies on numerical solutions. 
Before describing one method to numerically solve the 
equations, let us briefly comment on a powerful method 



to get an overview about the parameter values for which 
we might expect interesting behavior. 

Linear stability analysis 

This method starts by realizing that often, there is one 
rather simple solution to a continuum theory. Indeed, as 
is the case for the above equations, there is a spatially 
homogenous steady state, which is obtained by setting 
all derivatives equal to zero. This leaves us with an alge- 
braic set of equations, which can be solved numerically. 
The corresponding homogenous densities are denoted 
by cj, Co, n+, n^, and «c,o 

Your turn: Derive the equations for the homogenous 
steady state of Eqs. (30), (31), and (36)-(38). 

A real system is permanently subject to fluctuations. 
What happens to the spatially homogenous state if it is 
slightly perturbed? To study this question, we can write 
the distributions in the form c*{x) = Cq + Sc^{x) and analo- 
gously for c, n*, n, and «c- Inserting these expressions 
into the dynamic equations and keeping only terms in 
linear order in the perturbation, we arrive at a set of linear 
differential equations. 

Your turn: Derive the differential equations up to 
linear order in the perturbations (5c*, 5c , dn*, Sn, and 
dfii. What happens to the zeroth order terms? 

The linearized equations are easier to solve than the 
nonlinear dynamic equations, sometimes even analytically. 
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We do not want to do this explicitly here, as it requires 
some advanced mathematical techniques. Suffice it here to 
state that the linear stability analysis not only allows to 
identify interesting regions in parameter space, but often 
even yields some properties of the new state. In particular, 
it can indicate, whether the states that the system assumes 
under conditions of an unstable homogenous state is sta- 
tionary or oscillatory. However, only a nonlinear analysis 
can give the real answer. To this end one usually solves 
the dynamic equations numerically. 

Numerical solution 

There are many sophisticated methods to numerically 
solve non-linear partial differential equations. We will 
briefly introduce one very simple method that makes 
direct contact to the considerations in the mathematical 
interlude. While it may be slow, it works in a remark- 
ably large number of cases. 

In this numerical scheme, we discretize space and 
time - to some end we are thus going to reverse the 
continuum limit. Let us start by discretizing time. We 
only consider the density at discrete time points t = kAt, 
where At is the discretization length and k some integer. 
Then the derivative with respect to time is replaced by 



c((fe+ l)At) -c(feAt) 
At ■ 



Thus, given the time evolution equation 



-c=/(c(£),t). 



(39) 



(40) 



we can obtain c{{k + l)At) from c{kAt) by writing 

c((fe + 1) At) = c(feAt) + Atf (c(feAt), kAt). (41) 

In addition to time, the filament and nucleator densi- 
ties also depend on space. Let Ci(kAt) be the density at 
position iAx, where Ax is the spatial discretization 
length and / an integer with 0 <i < L/Ax and consider 
the equation 



— c[x, t) = ;(x, t). 



(42) 



An approximate expression for c,((A: + l)Ai) is then 
obtained from 

Cj((fe + l)At) = Cj{feAt) + At{ii{kAt) -ii^i{kAt))/Ax. (43) 

In this expression /, denotes the current from site i - 1 
to site i. Since CiAx is number of filaments in the bin i, 
this equation is the formal analog of Figure 2a. 

The discretized expressions for / = vc is 



if V > 0 and 



)i = VCi 



(45) 



if V < 0. Inserting this into Eq. (43), we get 

Ci((fe + 1) At) = Ci(feAt) + At{vci-i - VCi)/ Ax, (46) 

where we have considered the case v > 0. Since c, is a 
particle density, it must always be positive. From this 
condition, we get 1 - vAt/Ax > 0 as a sufficient condi- 
tion. That is, the values of At and Ax are not completely 
independent! Finally, we provide the expression for the 
diffusion current 



ji = D(ci_i - Ci)/Ax. 



(47) 



]i = VCi-i 



(44) 



Your turn: What is the condition on Ax and At 
resulting from this current? 

In Figure 5 we present space-time plots of the total 
NPF and filament densities in the case of an unstable 
homogenous distribution. In this case, the system indeed 
self-organizes into a traveling wave: Starting from a ran- 
dom initial condition, the filaments and NPFs form a 
distribution that moves at constant velocity. 

We can obtain additional insight into this state by 
plotting the various densities, see Figure 6. It turns out 
that there are practically only filaments of one orienta- 
tion, while the density of filaments of the opposite 
orientation is negligible. Similarly, this holds for the cor- 
responding filament-bound NPFs. Depending on the 
initial state, either one of the orientations will win and a 
wave either moving to the right or to the left will 
appear. Note, that we have not included any directional 
motion for the NPFs into our description. The apparent 
motion of the corresponding densities is a result of dif- 
fusion as well as binding to and unbinding from fila- 
ments: The peaks of the respective distributions of 
the NPFs and the filaments are shifted with the NPFs 
lagging behind, such that NPFs bind preferentially ahead 
of its maximum density. 

Conclusion 

In this tutorial, we have given an introduction into con- 
tinuum theories as a tool for analyzing cytoskeletal pro- 
cesses. We gave one example of a description that was 
derived from a discrete (microscopic) picture, and one 
example that was developed entirely within the frame of 
the continuum approach and was based on the continu- 
ity equation. Let us mention, that not all relevant quan- 
tities obey the continuity equation. If you consider, for 
example, a cytoskeletal network in two or three spatial 
dimensions, then one might need to take into account 
that, on large scales, the network can present orienta- 
tional order. Since the filament orientation is not a 
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Figure 5 Spontaneous cytoskeletal wave. Numerical solution to 
the dynamic equations (30)-(38). a) Filament density, b) Density of 
filament-bound NPFs. Warmer colors indicate higher densities. 
Parameter values are v = 0.1 ;jm/s, D = 0.01 /LjmVs, = 0.1 /Ljm^/s, 
v= 0.1 s"\ Vd = 0.1 s"\ cod = 0.1 s"\ coa = 0.01 /jm/s, and co, = 100 ^m^. 



conserved quantity, its dynamics is not given by the 
continuity equation. To derive dynamic equations for 
such fields one either relies on microscopic theories or 
on symmetry considerations. 

In the examples we discussed, we assumed a determi- 
nistic evolution of the system. In reality, the molecular 
processes underlying cellular processes are stochastic. In 
general, it is not an easy task to obtain the correct order 
of magnitude of stochastic effects. As long as one is only 
interested in qualitative behavior, it is often sufficient to 
capture molecular noise by effective diffusion terms. 

Probably a more serious point is that we neglected in 
our examples mechanical forces altogether. However, in 
cells and often also in reconstituted systems forces play an 
important role. These can be incorporated into continuum 
theories by accounting also for momentum conservation 
in addition to matter conservation. Sources and sinks are 




x(jwm) 

Figure 6 Filament and NPF densities in a wave Densities c*, n*, 
ric of plus-filaments (green), plus-NPFs (blue), and cytosolic NPFs 
(red) obtained from a numerical solution to the dynamic equations 
(30)-(38). Densities correspond to the state at time f = 100s in 
Figure 3. The densities of c and n' are too small (< lO"'') to be seen 
in the graph. 



in this case given by external forces, while the momentum 
flux density is given by the stress tensor. In addition to vis- 
cous or visco-elastic properties, the stress tensor also 
accounts for the stress generated by active processes [4]. 
In contrast, energy conservation is usually not an issue as 
cellular systems are embedded in a heat reservoir which 
maintains the system at constant temperature. 

In addition to the formal aspects of the descriptions 
developed above, let us comment on the biological sig- 
nificance of the processes considered. There is at least 
indirect experimental evidence for a feedback between 
the activity of nucleation promoting factors and actin 
filaments in neutrophils [14], although the molecular 
players involved are probably not completely known and 
mechanisms remain to be identified. As mentioned 
already above, the assembly dynamics of cytoskeletal 
filaments was also greatly simplified to not distract the 
reader from the focus of this tutorial. Indeed, neither 
the assembly nor the disassembly rates of filaments are 
expected to remain constant in time as typically they 
depend effectively on the filament length [19]. The 
assumption of constant assembly and disassembly rates 
prompted us to also assume a rate at which filaments 
spontaneously disintegrate. Although this is again unli- 
kely to be the case in living cells, the qualitative beha- 
vior of the artificial system is comparable to specific 
situations in more complicated scenarios [16]. A thor- 
ough discussion of realistic applications of continuum 
theories to cytoskeletal or other-cell-biological systems 
goes beyond the aim of the present tutorial. 

In summary, continuum descriptions provide powerful 
tools to analyze cellular processes. In a number of cases. 
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they have provided valuable information about biological 
and, in particular, cellular structures and processes. 
It will be exciting to see, whether they can also play a 
role in designing new therapeutic approaches. 
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